home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / ode-initval / test.c < prev   
Encoding:
C/C++ Source or Header  |  2002-04-18  |  15.0 KB  |  614 lines
  1. /* ode-initval/test_odeiv.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. /* Author:  G. Jungman
  21.  */
  22. #include <config.h>
  23. #include <stdlib.h>
  24. #include <stdio.h>
  25. #include <math.h>
  26. #include <gsl/gsl_test.h>
  27. #include <gsl/gsl_errno.h>
  28. #include <gsl/gsl_math.h>
  29. #include <gsl/gsl_matrix.h>
  30. #include <gsl/gsl_linalg.h>
  31. #include <gsl/gsl_ieee_utils.h>
  32. #include <gsl/gsl_odeiv.h>
  33.  
  34. int rhs_linear (double t, const double y[], double f[], void *params);
  35. int jac_linear (double t, const double y[], double *dfdy, double dfdt[],
  36.                 void *params);
  37. int rhs_sin (double t, const double y[], double f[], void *params);
  38. int jac_sin (double t, const double y[], double *dfdy, double dfdt[],
  39.              void *params);
  40. int rhs_exp (double t, const double y[], double f[], void *params);
  41. int jac_exp (double t, const double y[], double *dfdy, double dfdt[],
  42.              void *params);
  43. int rhs_stiff (double t, const double y[], double f[], void *params);
  44. int jac_stiff (double t, const double y[], double *dfdy, double dfdt[],
  45.                void *params);
  46. void test_stepper_linear (const gsl_odeiv_step_type * T, double h,
  47.                          double base_prec);
  48. void test_stepper_sin (const gsl_odeiv_step_type * T, double h, double base_prec);
  49. void test_stepper_exp (const gsl_odeiv_step_type * T, double h, double base_prec);
  50. void test_stepper_stiff (const gsl_odeiv_step_type * T, double h, double base_prec);
  51. void test_evolve_system_flat (gsl_odeiv_step * step,
  52.                               const gsl_odeiv_system * sys,
  53.                               double t0, double t1, double hstart,
  54.                               double y[], double yfin[],
  55.                               double err_target, const char *desc);
  56.  
  57. void test_evolve_system (const gsl_odeiv_step_type * T,
  58.                          const gsl_odeiv_system * sys,
  59.                          double t0, double t1, double hstart,
  60.                          double y[], double yfin[],
  61.                          double err_target, const char *desc);
  62.  
  63. void test_evolve_exp (const gsl_odeiv_step_type * T, double h, double err);
  64. void test_evolve_sin (const gsl_odeiv_step_type * T, double h, double err);
  65. void test_evolve_stiff1 (const gsl_odeiv_step_type * T, double h, double err);
  66. void test_evolve_stiff5 (const gsl_odeiv_step_type * T, double h, double err);
  67.  
  68. /* RHS for a + b t */
  69.  
  70. int
  71. rhs_linear (double t, const double y[], double f[], void *params)
  72. {
  73.   f[0] = 0.0;
  74.   f[1] = y[0];
  75.   return GSL_SUCCESS;
  76. }
  77.  
  78. int
  79. jac_linear (double t, const double y[], double *dfdy, double dfdt[],
  80.         void *params)
  81. {
  82.   gsl_matrix dfdy_mat;
  83.   dfdy_mat.size1 = 2;
  84.   dfdy_mat.size2 = 2;
  85.   dfdy_mat.tda = 2;
  86.   dfdy_mat.data = dfdy;
  87.   dfdy_mat.block = 0;
  88.   gsl_matrix_set (&dfdy_mat, 0, 0, 0.0);
  89.   gsl_matrix_set (&dfdy_mat, 0, 1, 0.0);
  90.   gsl_matrix_set (&dfdy_mat, 1, 0, 1.0);
  91.   gsl_matrix_set (&dfdy_mat, 1, 1, 0.0);
  92.   dfdt[0] = 0.0;
  93.   dfdt[1] = 0.0;
  94.   return GSL_SUCCESS;
  95. }
  96.  
  97. gsl_odeiv_system rhs_func_lin = {
  98.   rhs_linear,
  99.   jac_linear,
  100.   2,
  101.   0
  102. };
  103.  
  104.  
  105. /* RHS for sin(t),cos(t) */
  106.  
  107. int
  108. rhs_sin (double t, const double y[], double f[], void *params)
  109. {
  110.   f[0] = -y[1];
  111.   f[1] = y[0];
  112.   return GSL_SUCCESS;
  113. }
  114.  
  115. int
  116. jac_sin (double t, const double y[], double *dfdy, double dfdt[],
  117.      void *params)
  118. {
  119.   gsl_matrix dfdy_mat;
  120.   dfdy_mat.data = dfdy;
  121.   dfdy_mat.size1 = 2;
  122.   dfdy_mat.size2 = 2;
  123.   dfdy_mat.tda = 2;
  124.   dfdy_mat.block = 0;
  125.   gsl_matrix_set (&dfdy_mat, 0, 0, 0.0);
  126.   gsl_matrix_set (&dfdy_mat, 0, 1, -1.0);
  127.   gsl_matrix_set (&dfdy_mat, 1, 0, 1.0);
  128.   gsl_matrix_set (&dfdy_mat, 1, 1, 0.0);
  129.   dfdt[0] = 0.0;
  130.   dfdt[1] = 0.0;
  131.   return GSL_SUCCESS;
  132. }
  133.  
  134. gsl_odeiv_system rhs_func_sin = {
  135.   rhs_sin,
  136.   jac_sin,
  137.   2,
  138.   0
  139. };
  140.  
  141.  
  142. /* RHS for a exp(t)+ b exp(-t) */
  143.  
  144. int
  145. rhs_exp (double t, const double y[], double f[], void *params)
  146. {
  147.   f[0] = y[1];
  148.   f[1] = y[0];
  149.   return GSL_SUCCESS;
  150. }
  151.  
  152. int
  153. jac_exp (double t, const double y[], double *dfdy, double dfdt[],
  154.      void *params)
  155. {
  156.   gsl_matrix dfdy_mat;
  157.   dfdy_mat.data = dfdy;
  158.   dfdy_mat.size1 = 2;
  159.   dfdy_mat.size2 = 2;
  160.   dfdy_mat.tda = 2;
  161.   dfdy_mat.block = 0;
  162.   gsl_matrix_set (&dfdy_mat, 0, 0, 0.0);
  163.   gsl_matrix_set (&dfdy_mat, 0, 1, 1.0);
  164.   gsl_matrix_set (&dfdy_mat, 1, 0, 1.0);
  165.   gsl_matrix_set (&dfdy_mat, 1, 1, 0.0);
  166.   dfdt[0] = 0.0;
  167.   dfdt[1] = 0.0;
  168.   return GSL_SUCCESS;
  169. }
  170.  
  171. gsl_odeiv_system rhs_func_exp = {
  172.   rhs_exp,
  173.   jac_exp,
  174.   2,
  175.   0
  176. };
  177.  
  178.  
  179. /* RHS for stiff example */
  180.  
  181. int
  182. rhs_stiff (double t, const double y[], double f[], void *params)
  183. {
  184.   f[0] = 998.0 * y[0] + 1998.0 * y[1];
  185.   f[1] = -999.0 * y[0] - 1999.0 * y[1];
  186.   return GSL_SUCCESS;
  187. }
  188.  
  189. int
  190. jac_stiff (double t, const double y[], double *dfdy, double dfdt[],
  191.        void *params)
  192. {
  193.   gsl_matrix dfdy_mat;
  194.   dfdy_mat.data = dfdy;
  195.   dfdy_mat.size1 = 2;
  196.   dfdy_mat.size2 = 2;
  197.   dfdy_mat.tda = 2;
  198.   dfdy_mat.block = 0;
  199.   gsl_matrix_set (&dfdy_mat, 0, 0, 998.0);
  200.   gsl_matrix_set (&dfdy_mat, 0, 1, 1998.0);
  201.   gsl_matrix_set (&dfdy_mat, 1, 0, -999.0);
  202.   gsl_matrix_set (&dfdy_mat, 1, 1, -1999.0);
  203.   dfdt[0] = 0.0;
  204.   dfdt[1] = 0.0;
  205.   return GSL_SUCCESS;
  206. }
  207.  
  208. gsl_odeiv_system rhs_func_stiff = {
  209.   rhs_stiff,
  210.   jac_stiff,
  211.   2,
  212.   0
  213. };
  214.  
  215.  
  216. void
  217. test_stepper_linear (const gsl_odeiv_step_type * T, double h,
  218.              double base_prec)
  219. {
  220.   int s = 0;
  221.   double y[2];
  222.   double yerr[2];
  223.   double t;
  224.   double del;
  225.   double delmax = 0.0;
  226.   int count = 0;
  227.  
  228.   gsl_odeiv_step *stepper = gsl_odeiv_step_alloc (T, 2);
  229.  
  230.   y[0] = 1.0;
  231.   y[1] = 0.0;
  232.  
  233.   for (t = 0.0; t < 4.0; t += h)
  234.     {
  235.       gsl_odeiv_step_apply (stepper, t, h, y, yerr, 0, 0, &rhs_func_lin);
  236.       del = fabs ((y[1] - (t + h)) / y[1]);
  237.       delmax = GSL_MAX_DBL (del, delmax);
  238.       if (del > (count + 1.0) * base_prec)
  239.     {
  240.       printf ("  LINEAR(%20.17g)  %20.17g  %20.17g  %8.4g\n", t + h, y[1],
  241.           t + h, del);
  242.       s++;
  243.     }
  244.       count++;
  245.     }
  246.  
  247.   gsl_test (s, "%s, linear [0,4], max relative error = %g",
  248.         gsl_odeiv_step_name (stepper), delmax);
  249.  
  250.   gsl_odeiv_step_free (stepper);
  251. }
  252.  
  253.  
  254. void
  255. test_stepper_sin (const gsl_odeiv_step_type * T, double h, double base_prec)
  256. {
  257.   int s = 0;
  258.   double y[2];
  259.   double yerr[2];
  260.   double t;
  261.   double del;
  262.   double delmax = 0.0;
  263.   int count = 0;
  264.  
  265.   gsl_odeiv_step *stepper = gsl_odeiv_step_alloc (T, 2);
  266.  
  267.   y[0] = 1.0;
  268.   y[1] = 0.0;
  269.  
  270.   for (t = 0.0; t < M_PI; t += h)
  271.     {
  272.       int stat;
  273.       double sin_th = sin (t + h);
  274.       gsl_odeiv_step_apply (stepper, t, h, y, yerr, 0, 0, &rhs_func_sin);
  275.       del = fabs ((y[1] - sin_th) / sin_th);
  276.       delmax = GSL_MAX_DBL (del, delmax);
  277.       {
  278.     if (t < 0.5 * M_PI)
  279.       {
  280.         stat = (del > (count + 1.0) * base_prec);
  281.       }
  282.     else if (t < 0.7 * M_PI)
  283.       {
  284.         stat = (del > 1.0e+04 * base_prec);
  285.       }
  286.     else if (t < 0.9 * M_PI)
  287.       {
  288.         stat = (del > 1.0e+06 * base_prec);
  289.       }
  290.     else
  291.       {
  292.         stat = (del > 1.0e+09 * base_prec);
  293.       }
  294.     if (stat != 0)
  295.       {
  296.         printf ("  SIN(%22.18g)  %22.18g  %22.18g  %10.6g\n", t + h, y[1],
  297.             sin_th, del);
  298.       }
  299.     s += stat;
  300.       }
  301.       count++;
  302.     }
  303.   if (delmax > 1.0e+09 * base_prec)
  304.     {
  305.       s++;
  306.       printf ("  SIN(0 .. M_PI)  delmax = %g\n", delmax);
  307.     }
  308.  
  309.   gsl_test (s, "%s, sine [0,pi], max relative error = %g",
  310.         gsl_odeiv_step_name (stepper), delmax);
  311.  
  312.   delmax = 0.0;
  313.   for (; t < 100.5 * M_PI; t += h)
  314.     {
  315.       gsl_odeiv_step_apply (stepper, t, h, y, yerr, 0, 0, &rhs_func_sin);
  316.       del = fabs (y[1] - sin (t));
  317.       delmax = GSL_MAX_DBL (del, delmax);
  318.       count++;
  319.     }
  320.  
  321.   if (del > count * 2.0 * base_prec)
  322.     {
  323.       s++;
  324.       printf ("  SIN(%22.18g)  %22.18g  %22.18g  %10.6g\n", t + h, y[1],
  325.           sin (t), del);
  326.     }
  327.  
  328.   gsl_test (s, "%s, sine [pi,100.5*pi], max absolute error = %g",
  329.         gsl_odeiv_step_name (stepper), delmax);
  330.  
  331.   gsl_odeiv_step_free (stepper);
  332. }
  333.  
  334.  
  335. void
  336. test_stepper_exp (const gsl_odeiv_step_type * T, double h, double base_prec)
  337. {
  338.   int s = 0;
  339.   double y[2];
  340.   double yerr[2];
  341.   double t;
  342.   double del, delmax = 0.0;
  343.   int count = 0;
  344.  
  345.   gsl_odeiv_step *stepper = gsl_odeiv_step_alloc (T, 2);
  346.  
  347.   y[0] = 1.0;
  348.   y[1] = 1.0;
  349.  
  350.   for (t = 0.0; t < 20.0; t += h)
  351.     {
  352.       double ex = exp (t + h);
  353.       gsl_odeiv_step_apply (stepper, t, h, y, yerr, 0, 0, &rhs_func_exp);
  354.       del = fabs ((y[1] - ex) / y[1]);
  355.       delmax = GSL_MAX_DBL (del, delmax);
  356.       if (del > (count + 1.0) * 2.0 * base_prec)
  357.     {
  358.       printf ("  EXP(%20.17g)  %20.17g  %20.17g  %8.4g\n", t + h, y[1],
  359.           ex, del);
  360.       s++;
  361.     }
  362.       count++;
  363.     }
  364.  
  365.   gsl_test (s, "%s, exponential [0,20], max relative error = %g",
  366.         gsl_odeiv_step_name (stepper), delmax);
  367.  
  368.   gsl_odeiv_step_free (stepper);
  369. }
  370.  
  371. void
  372. test_stepper_stiff (const gsl_odeiv_step_type * T, double h, double base_prec)
  373. {
  374.   int s = 0;
  375.   double y[2];
  376.   double yerr[2];
  377.   double t;
  378.   double del, delmax = 0.0;
  379.   int count = 0;
  380.  
  381.   gsl_odeiv_step *stepper = gsl_odeiv_step_alloc (T, 2);
  382.  
  383.   y[0] = 1.0;
  384.   y[1] = 0.0;
  385.  
  386.   for (t = 0.0; t < 20.0; t += h)
  387.     {
  388.       gsl_odeiv_step_apply (stepper, t, h, y, yerr, NULL, NULL,
  389.                 &rhs_func_stiff);
  390.  
  391.       if (t > 0.04)
  392.     {
  393.       double arg = t + h;
  394.       double e1 = exp (-arg);
  395.       double e2 = exp (-1000.0 * arg);
  396.       double u = 2.0 * e1 - e2;
  397.       /* double v = -e1 + e2; */
  398.       del = fabs ((y[0] - u) / y[0]);
  399.       delmax = GSL_MAX_DBL (del, delmax);
  400.  
  401.       if (del > (count + 1.0) * 100.0 * base_prec)
  402.         {
  403.           printf ("  STIFF(%20.17g)  %20.17g  %20.17g  %8.4g\n", arg,
  404.               y[0], u, del);
  405.           s++;
  406.         }
  407.     }
  408.       count++;
  409.     }
  410.  
  411.   gsl_test (s, "%s, stiff [0,20], max relative error = %g",
  412.         gsl_odeiv_step_name (stepper), delmax);
  413.  
  414.   gsl_odeiv_step_free (stepper);
  415. }
  416.  
  417. void
  418. test_evolve_system_flat (gsl_odeiv_step * step,
  419.              const gsl_odeiv_system * sys,
  420.              double t0, double t1, double hstart,
  421.              double y[], double yfin[],
  422.              double err_target, const char *desc)
  423. {
  424.   int s = 0;
  425.   double frac;
  426.  
  427.   double t = t0;
  428.   double h = hstart;
  429.  
  430.   gsl_odeiv_evolve *e = gsl_odeiv_evolve_alloc (sys->dimension);
  431.  
  432.   while (t < t1)
  433.     {
  434.       gsl_odeiv_evolve_apply (e, NULL, step, sys, &t, t1, &h, y);
  435.     }
  436.  
  437.   frac = fabs ((y[1] - yfin[1]) / yfin[1]) + fabs ((y[0] - yfin[0]) / yfin[0]);
  438.  
  439.   if (frac > 2.0 * e->count * err_target)
  440.     {
  441.       printf ("FLAT t = %.5e  y0 = %g y1= %g\n", t, y[0], y[1]);
  442.       s++;
  443.     }
  444.  
  445.   gsl_test (s, "%s, %s, evolve, no control, max relative error = %g",
  446.         gsl_odeiv_step_name (step), desc, frac);
  447.  
  448.   gsl_odeiv_evolve_free (e);
  449. }
  450.  
  451.  
  452. void
  453. test_evolve_system (const gsl_odeiv_step_type * T,
  454.             const gsl_odeiv_system * sys,
  455.             double t0, double t1, double hstart,
  456.             double y[], double yfin[],
  457.             double err_target, const char *desc)
  458. {
  459.   int s = 0;
  460.   double frac;
  461.  
  462.   double t = t0;
  463.   double h = hstart;
  464.  
  465.   gsl_odeiv_step * step = gsl_odeiv_step_alloc (T, sys->dimension);
  466.  
  467.   gsl_odeiv_control *c =
  468.     gsl_odeiv_control_standard_new (0.0, err_target, 1.0, 1.0);
  469.   gsl_odeiv_evolve *e = gsl_odeiv_evolve_alloc (sys->dimension);
  470.  
  471.   while (t < t1)
  472.     {
  473.       gsl_odeiv_evolve_apply (e, c, step, sys, &t, t1, &h, y);
  474.       /* printf ("SYS  t = %.18e h = %g y0 = %g y1= %g\n", t, h, y[0], y[1]); */
  475.     }
  476.  
  477.   frac = fabs ((y[1] - yfin[1]) / yfin[1]) + fabs ((y[0] - yfin[0]) / yfin[0]);
  478.  
  479.   if (frac > 2.0 * e->count * err_target)
  480.     {
  481.       printf ("SYS  t = %.5e h = %g y0 = %g y1= %g\n", t, h, y[0], y[1]);
  482.       s++;
  483.     }
  484.  
  485.   gsl_test (s, "%s, %s, evolve, standard control, relative error = %g",
  486.         gsl_odeiv_step_name (step), desc, frac);
  487.  
  488.   gsl_odeiv_evolve_free (e);
  489.   gsl_odeiv_control_free (c);
  490.   gsl_odeiv_step_free (step);
  491. }
  492.  
  493.  
  494. void
  495. test_evolve_exp (const gsl_odeiv_step_type * T, double h, double err)
  496. {
  497.   double y[2];
  498.   double yfin[2];
  499.  
  500.   y[0] = 1.0;
  501.   y[1] = 1.0;
  502.   yfin[0] = exp (10.0);
  503.   yfin[1] = yfin[0];
  504.   test_evolve_system (T, &rhs_func_exp, 0.0, 10.0, h, y, yfin, err, "exp [0,10]");
  505. }
  506.  
  507. void
  508. test_evolve_sin (const gsl_odeiv_step_type * T, double h, double err)
  509. {
  510.   double y[2];
  511.   double yfin[2];
  512.  
  513.   y[0] = 1.0;
  514.   y[1] = 0.0;
  515.   yfin[0] = cos (2.0);
  516.   yfin[1] = sin (2.0);
  517.   test_evolve_system (T, &rhs_func_sin, 0.0, 2.0, h, y, yfin, err, "sine [0,2]");
  518. }
  519.  
  520. void
  521. test_evolve_stiff1 (const gsl_odeiv_step_type * T, double h, double err)
  522. {
  523.   double y[2];
  524.   double yfin[2];
  525.  
  526.   y[0] = 1.0;
  527.   y[1] = 0.0;
  528.   {
  529.     double arg = 1.0;
  530.     double e1 = exp (-arg);
  531.     double e2 = exp (-1000.0 * arg);
  532.     yfin[0] = 2.0 * e1 - e2;
  533.     yfin[1] = -e1 + e2;
  534.   }
  535.   test_evolve_system (T, &rhs_func_stiff, 0.0, 1.0, h, y, yfin, err,
  536.                       "stiff [0,1]");
  537. }
  538.  
  539. void
  540. test_evolve_stiff5 (const gsl_odeiv_step_type * T, double h, double err)
  541. {
  542.   double y[2];
  543.   double yfin[2];
  544.  
  545.   y[0] = 1.0;
  546.   y[1] = 0.0;
  547.   {
  548.     double arg = 5.0;
  549.     double e1 = exp (-arg);
  550.     double e2 = exp (-1000.0 * arg);
  551.     yfin[0] = 2.0 * e1 - e2;
  552.     yfin[1] = -e1 + e2;
  553.   }
  554.   test_evolve_system (T, &rhs_func_stiff, 0.0, 5.0, h, y, yfin, err,
  555.                       "stiff [0,5]");
  556. }
  557.  
  558.  
  559.  
  560. int
  561. main (void)
  562. {
  563.   int i;
  564.  
  565.   struct ptype
  566.   {
  567.     const gsl_odeiv_step_type *type;
  568.     double h;
  569.   }
  570.   p[20];
  571.  
  572.   p[0].type = gsl_odeiv_step_rk2;
  573.   p[0].h = 1.0e-03;
  574.   p[1].type = gsl_odeiv_step_rk2imp;
  575.   p[1].h = 1.0e-03;
  576.   p[2].type = gsl_odeiv_step_rk4;
  577.   p[2].h = 1.0e-03;
  578.   p[3].type = gsl_odeiv_step_rk4imp;
  579.   p[3].h = 1.0e-03;
  580.   p[4].type = gsl_odeiv_step_rkf45;
  581.   p[4].h = 1.0e-03;
  582.   p[5].type = gsl_odeiv_step_rk8pd;
  583.   p[5].h = 1.0e-03;
  584.   p[6].type = gsl_odeiv_step_rkck;
  585.   p[6].h = 1.0e-03;
  586.   p[7].type = gsl_odeiv_step_bsimp;
  587.   p[7].h = 0.1;
  588.   p[8].type = gsl_odeiv_step_gear1;
  589.   p[8].h = 1.0e-03;
  590.   p[9].type = gsl_odeiv_step_gear2;
  591.   p[9].h = 1.0e-03;
  592.   p[10].type = 0;
  593.  
  594.   gsl_ieee_env_setup ();
  595.  
  596.   for (i = 0; p[i].type != 0; i++)
  597.     {
  598.       test_stepper_linear (p[i].type, p[i].h, GSL_SQRT_DBL_EPSILON);
  599.       test_stepper_sin (p[i].type, p[i].h / 10.0, GSL_SQRT_DBL_EPSILON);
  600.       test_stepper_exp (p[i].type, p[i].h / 10.0, GSL_SQRT_DBL_EPSILON);
  601.       test_stepper_stiff (p[i].type, p[i].h / 10.0, GSL_SQRT_DBL_EPSILON);
  602.     }
  603.  
  604.   for (i = 0; p[i].type != 0; i++)
  605.     {
  606.       test_evolve_exp (p[i].type, p[i].h, GSL_SQRT_DBL_EPSILON);
  607.       test_evolve_sin (p[i].type, p[i].h, GSL_SQRT_DBL_EPSILON);
  608.       test_evolve_stiff1 (p[i].type, p[i].h, 1e-5);
  609.       test_evolve_stiff5 (p[i].type, p[i].h, 1e-5);
  610.     }
  611.  
  612.   exit (gsl_test_summary ());
  613. }
  614.